/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 *
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

package models.com.javable.cese.models;

import java.io.Serializable;

import com.javable.cese.templates.AbstractAPModel;

/**
 * The FitzHugh-Nagumo Simplified Cardiac Action Potential Model, 1961. Detailed
 * list of equations and model description are provided in
 * <p>
 * Biophys J 1961;1:445-466,<br>
 * Chaos, Solitons and Fractals 1996;7:293-301.
 */

public class FitzHughNagumoModel_1961 extends AbstractAPModel implements
        Serializable {

    /*
     * The beginning here lists all the variables and defines them (this code
     * uses all global variables)
     */

    protected double u = 0.0; // Dimensionless potential

    public void setu(double newu) {
        u = newu;
    }

    public double getu() {
        return u;
    }

    protected double vr = -85.0; // Resting potential (mV)

    public void setvr(double newvr) {
        vr = newvr;
    }

    public double getvr() {
        return vr;
    }

    protected double vp = 15.0; // Peak potential (mV)

    public void setvp(double newvp) {
        vp = newvp;
    }

    public double getvp() {
        return vp;
    }

    protected double alpha = 0.0; // Excitation threshold

    public void setalpha(double newalpha) {
        alpha = newalpha;
    }

    public double getalpha() {
        return alpha;
    }

    protected double vth = -75.0; // Excitation threshold potential (mV)

    public void setvth(double newvth) {
        vth = newvth;
    }

    public double getvth() {
        return vth;
    }

    protected double c1 = 0.175; // Time constant 1 (ms^-1)

    public void setc1(double newc1) {
        c1 = newc1;
    }

    public double getc1() {
        return c1;
    }

    protected double c2 = 0.03; // Time constant 2 (ms^-1)

    public void setc2(double newc2) {
        c2 = newc2;
    }

    public double getc2() {
        return c2;
    }

    protected double vv = 0.0; // Slow current

    public void setvv(double newvv) {
        vv = newvv;
    }

    public double getvv() {
        return vv;
    }

    protected double b = 0.011; // Time constant of the slow current (ms^-1)

    public void setb(double newb) {
        b = newb;
    }

    public double getb() {
        return b;
    }

    protected double d = 0.55; // Slow current resistance

    public void setd(double newd) {
        d = newd;
    }

    public double getd() {
        return d;
    }

    protected double cm = 0.01; // Cell capacitance (pF)

    public void setcm(double newcm) {
        cm = newcm;
    }

    public double getcm() {
        return cm;
    }

    private double iont; // Total Ion Flow (mM/ms)

    public double getiont() {
        return iont;
    }

    /* ** CONSTRUCTOR ** */
    public FitzHughNagumoModel_1961() {
        super();

        setst(-1.0);
    }

    public void setCurrents() {
        component_membrane();
        component_recovery_variable();

        currentIt();
    }

    public void doSimulate() {
        setCurrents();

        vnew = v - (it / cm) * dt;
        dvdtnew = (vnew - v) / dt;
        t = t + dt;
        v = vnew;
        dvdt = dvdtnew;
    }

    /** ***************************************************** */

    /* Functions that describe the currents begin here */

    public void component_membrane() {
        u = (v - vr) / (vp - vr);
    }

    public void component_recovery_variable() {
        vv = vv + (dt * (b * (u - (d * vv))));
    }

    /*
     * Total sum of currents is calculated here, if the time is between stimtime =
     * 0 and stimtime = tstim, a stimulus is applied
     */

    public void currentIt() // Calculates Total Current
    {
        alpha = (vth - vr) / (vp - vr);
        iont = ((c1 * u * (u - alpha) * (u - 1.0)) + (c2 * vv));

        if ((t >= tstim) && (t < (tstim + dt))) {
            stimtime = 0;
            tstim = tstim + bcl;
            dvdtmax = false;
        }

        if ((stimtime >= 0) && (stimtime < stl)) {
            it = st + iont;
        } else {
            it = iont;
        }

        stimtime = stimtime + dt;
    }

}